Evaluate and detrend OISST anomaly by 30yrs mean

R code modified from https://github.com/mjacox/Thermal_Displacement/blob/master/oisst_an.m

First calculate 1982-2011 SST monthly mean. To evaluate SST anomaly:

1. SST.1 = Monthly SST - SST_30yr_mean

2. SST.anomaly = SST.1 - Linear trend by monthly series from 1982-202205

3. Land mask from R/01-2_OISST_land_seaice_mask.Rmd

## The same code from R/01-2_OISST_land_seaice_mask.Rmd(and can check that testing figure)
land_mask <- fread("unzip -p ../data_src/mapfiles/sea_icemask_025d.csv.zip")[, .(lon, lat, x, y, landmask)]
setorder(land_mask, -y, x)
land_mask[,ry:=rowid(x)]
yy <- unique(land_mask[,.(y, ry)]) ## check, NOTE that in OISST.nc file, y is from Northest (89.875) to Southest (-89.875)

landm <- dcast(land_mask, x ~ ry, value.var = "landmask")
xx <- as.numeric(landm[,1]$x)
landm <- as.matrix(landm[,-1])

clim_years = seq(1982,2011) #for climatology
climyrs<- clim_years[length(clim_years)] - clim_years[1] + 1 #30 yrs
monstr <- c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec")

Recalculate_30yrs = FALSE

if (Recalculate_30yrs) {
  plan(multisession)
  options(future.globals.maxSize= 1048576000)

  stm <- future_lapply(1:12, function(j) {
    monj <- fifelse(j<10, paste0("0",j), paste0(j))
    jstr <- monstr[j]
    for (i in clim_years) {
      x <- read_stars(paste0("../data_src/oisst/monthly_sst/", i, monj, "_sst.nc"))
      names(x)[1] <- jstr 
      x[[1]][which(landm==1)] <- NA_real_
      ice <- read_stars(paste0("../data_src/oisst/monthly_icemask/", i, monj, "_icemask.nc"))
      x[[1]][which(ice[[1]]==1)] <- NA_real_
      if (i == clim_years[1]) {
        stmx <- x
        datex<- as.Date(paste0(i, monj, "01"), format="%Y%m%d")
      } else {
        stmx <- c(stmx, x)
        datex<- c(datex, as.Date(paste0(i, monj, "01"), format="%Y%m%d"))
      }
    }
    names(stmx) <- rep(jstr, length(names(stmx)))
    stmx <- merge(stmx) %>% st_set_dimensions(3, values = as.POSIXct(datex), names = "time") %>% 
        aggregate(by=paste0(climyrs, " years"), FUN=mean, na.rm=TRUE)
    names(stmx)[1] <- jstr
    stmx <- stmx %>% select(jstr) %>% adrop
    return (stmx)
  })
}

if (Recalculate_30yrs) {
  # Note that in an earlier version used "GLDASp4 land mask" but current version use "GLDASp5"
  # The differnt version of land mask cause difference in the definition of land areas, so that sea-mask
  # and then also cause difference in climatology mean result
  save(stm, file="../data_src/stats/sst_1982_2011month_climatology_nc.RData")
} else {
  load("../data_src/stats/sst_1982_2011month_climatology_nc.RData")
}

stm

library(grid)

plotx <- vector("list", length = 12)
pstrx <- '';
for (j in 1:12) {
  plotx[[j]] <- gplotx(stm[[j]], monstr[j], returnx = TRUE, minz = -2.5, maxz = 34.5)
  pstrx <- paste0(pstrx, "plotx[[", j, "]],") 
}

lay1 <- rbind(c(1,1,2,2,3,3),
              c(4,4,5,5,6,6),
              c(7,7,8,8,9,9),
              c(10,10,11,11,12,12))
evplot <- paste0('grid.arrange(', pstrx, ' layout_matrix=lay1)')

gx <- eval(parse(text=evplot))

To double check if the plots of climatology are right, use existed data that download from NOAA (NetCDF OISST, 0.5d, monthly mean)

• Metadata: https://psl.noaa.gov/data/gridded/data.noaa.oisst.v2.html

• https://downloads.psl.noaa.gov/Datasets/noaa.oisst.v2/sst.mnmean.nc (downloaded 202206 and the data is 1981/12 to 2022/05)

if (!require("ncdf4")) install.packages("ncdf4"); library(ncdf4)

sst_mnfile = "../data_src/stats/sst.mnmean.nc"
if (!file.exists(sst_mnfile)) {
  tryCatch({
    curl::curl_download("https://downloads.psl.noaa.gov/Datasets/noaa.oisst.v2/sst.mnmean.nc", destfile = sst_mnfile)
  }, error = function(e) paste0(e))
} 

nx0 <- nc_open(sst_mnfile) 
print(nx0) ## sst[lon,lat,time]   (Chunking: [360,180,464])
latn1<- ncvar_get(nx0, "lat") # 89.5 - -89.5
lngn1<- ncvar_get(nx0, "lon") # 0 - 359.5
time<- ncvar_get(nx0, "time") # 1981-12.01 - 2020-07-01, we now check every March, August, and December
date<- time %>%  as.Date(origin="1800-01-01 00:00:0.0") 

maridx <- seq(4, 464, by = 12) #check date[maridx]
augidx <- seq(9, 464, by = 12)
decidx <- seq(13, 464, by = 12)

get_sstmx <- function (midx, minz=-2.5, maxz=34.5) {
  dtx <- as.data.table(ncvar_get(nx0, "sst")[,,midx]) %>% 
    .[,.(sstm=mean(value, na.rm=TRUE)), by=.(V1,V2)] %>%
    .[,`:=`(longitude=fifelse(lngn1[V1]>180, lngn1[V1]-360, lngn1[V1]), latitude=latn1[V2])] %>%
    .[,.(longitude, latitude, sstm)]
  setnames(dtx, 1:2 , c("x", "y"))
  
  gx <- gplotx(dtx, var="sstm", returnx=TRUE, minz=minz, maxz=maxz, no_coord=TRUE, maxlim=180) 
  gx <- gx +
  #gx <- ggplot() + 
  #  geom_tile(data=dt, aes_string(x="longitude", y="latitude", fill="sstm"), alpha=0.8) + 
  #  scale_fill_viridis(limits=c(minsst, maxsst)) +
     geom_sf(data = ne_coastline(scale = "large", returnclass = "sf"), color = 'darkgray', size = .3) #+
  #  coord_sf() + 
  #  xlim(c(-180, 180)) + ylim(c(-90, 90))
  
  return(gx)  
}

gx3 <- get_sstmx(maridx)
gx8 <- get_sstmx(augidx)
gx12 <-get_sstmx(decidx)

gt3 <- gplotx(as.data.table(stm[[3]]) %>% .[,`:=`(x=fifelse(x>180, x-360, x))], "Mar", returnx=TRUE, minz=-2.5, maxz=34.5, maxlim=180)  
gt8 <- gplotx(as.data.table(stm[[8]]) %>% .[,`:=`(x=fifelse(x>180, x-360, x))], "Aug", returnx=TRUE, minz=-2.5, maxz=34.5, maxlim=180)  
gt12 <- gplotx(as.data.table(stm[[12]]) %>% .[,`:=`(x=fifelse(x>180, x-360, x))], "Dec", returnx=TRUE, minz=-2.5, maxz=34.5, maxlim=180)  

lay2 <- rbind(c(1,1,2,2),
              c(3,3,4,4),
              c(5,5,6,6))

grid.arrange(gx3, gt3, gx8, gt8, gx12, gt12, layout_matrix=lay2)

i <- 2022
j <- 4
monj <- fifelse(j<10, paste0("0",j), paste0(j))
jstr <- monstr[j]
  
x <- read_stars(paste0("../data_src/oisst/monthly_sst/", i, monj, "_sst.nc"))
names(x)[1] <- "anom" 
#dim(x[[1]]) #1440. 720
#dim(landm)  #1440, 720
#dim(ice)    #1440, 720
x[[1]][which(landm==1)] <- NA_real_
ice <- read_stars(paste0("../data_src/oisst/monthly_icemask/", i, monj, "_icemask.nc"))
x[[1]][which(ice[[1]]==1)] <- NA_real_

#xt <- matrix()
#length(xt) <-1440 * 720
#dim(xt) <- c(1440, 720)
## if one of x and stm is NA, then the substraction is NA (so numbers of NA values increases)
x[[1]] <- x[[1]] - stm[[j]][[1]]

## Just check if plot 2022-04
gx <- gplotx(x, "anom", minz=-3.5, maxz=3.5, return=TRUE) #This anomaly is substrat long-term mean (30yrs, 1982-2011)
y <- read_stars("../data_src/oisst/monthly_anom/202204_anom.nc")
names(y)[1] <- "anom" 
gy <- gplotx(y, "anom", minz=-3.5, maxz=3.5, return=TRUE)

layt <- rbind(c(1,1,2,2))
grid.arrange(gx, gy, layout_matrix=layt)

yrng <- seq(1982,2022)
endt <- "2022-05-22"
trackdate <- seq.Date(as.IDate(endt)-6, as.IDate(endt), by="day")
curryr <- year(as.IDate(endt))
currmo <- month(as.IDate(endt))

j=4 #just testing, arbitrary select a month
monj <- fifelse(j<10, paste0("0",j), paste0(j))
jstr <- monstr[j]
for (i in yrng) {
  if (i==curryr & j>(currmo-1)) break
  z <- read_stars(paste0("../data_src/oisst/monthly_anom_icemask/", i, monj, "_anom.nc"))
  if (i==yrng[1]) {
    stdrx <- z
    datex<- as.Date(paste0(i, monj, "01"), format="%Y%m%d")
  } else {
    stdrx <- c(stdrx, z)
    datex<- c(datex, as.Date(paste0(i, monj, "01"), format="%Y%m%d"))
  }
}
names(stdrx) <- rep(jstr, length(names(stdrx)))
xt <- merge(stdrx)
st_set_dimensions(xt, 3, values = as.POSIXct(datex), names = "time")

## stars object with 3 dimensions and 1 attribute
## attribute(s), summary of first 1e+05 cells:
##       Min.    1st Qu.      Median       Mean    3rd Qu.     Max.  NA's
## X  -2.8777 -0.4862444 -0.09030438 -0.2813758 0.05061904 1.245278 95907
## dimension(s):
##      from   to offset delta  refsys point
## x       1 1440      0  0.25      NA    NA
## y       1  720     90 -0.25      NA    NA
## time    1   41     NA    NA POSIXct    NA
##                                           values x/y
## x                                           NULL [x]
## y                                           NULL [y]
## time 1982-04-01 08:00:00,...,2022-04-01 08:00:00

seqx <- rbind(c(`377`= as.numeric(xt[[1]][377,41,])), ## arbitrary select a point, check which(!is.na(xt[[1]][,41,]))
              c(`1112`= as.numeric(xt[[1]][1112,41,])),
              c(`test`= sort(rnorm(41)) + rnorm(41)))
colnames(seqx) <- gsub("-", "_", datex)
rownames(seqx) = c("0377_41", "1222_41", "test")
yhat <- stats::predict(lm(seqx[3,] ~ seq_along(datex)))
seqx[1, 28:30] <- NA_real_ #insert NA
seqx[2, 32:34] <- NA_real_ #insert NA
predy <- function(x) { return (stats::predict(lm(x ~ seq_along(datex)))) }
haty <- apply(seqx, MARGIN=1, FUN=function(y){ #predy use predict will auto exclude NA, so that out length is shorter
    if (all(is.na(y))) return(y)
    yh <- as.numeric(predy(y)) 
    y[!is.na(y)] <- yh
    return(c(as.numeric(y)))
}) %>% as.data.frame()  
all.equal(haty$test, as.numeric(yhat)) #TRUE

## [1] TRUE

#test using package pracma detrend function
if (!require("pracma")) install.packages("pracma"); library(pracma)

## Loading required package: pracma

## 
## Attaching package: 'pracma'

## The following objects are masked from 'package:magrittr':
## 
##     and, mod, or

ydet <- as.numeric(detrend(seqx[3,], 'linear'))
all.equal(ydet, as.numeric(seqx[3,]) - as.numeric(yhat)) #TRUE

## [1] TRUE

#test using package pracma detrend function
if (!require("pracma")) install.packages("pracma"); library(pracma)

plot_pred <- function(x, y, yhat=NULL) {
  if (is.null(yhat)) {
    yhat <- stats::predict(lm(y ~ x))
  } else {
    yh <- stats::predict(lm(y ~ x))
    if (!isTRUE(all.equal(as.numeric(yh), yhat))) {
      print(paste0("Not equal yhat with lm: ", paste(yh, collapse=",")))
    }
  }
  plot(x=x, y=y)
  abline(lm(y ~ x))
  points(x=x, y=yhat, col="red", pch=4)
  points(x=x, y=c(y-yhat), col="green", pch=13)
}

plot_pred(seq_along(datex), seqx[3,], haty$`test`)

save(seqx, file="../data_src/stats/test_stars_predict00.RData")

Recalculate_anom <- FALSE
if (Recalculate_anom) {
  
for (i in yrng) {
  mmx <- fifelse(i==curryr, currmo-1L, 12L)
  for (j in 1:mmx) {
    monj <- fifelse(j<10, paste0("0",j), paste0(j))
    jstr <- monstr[j]
    filet <- paste0("../data_src/oisst/monthly_anom_icemask/", i, monj, "_anom.nc")
    if (!file.exists(filet)) {
      #z<- read_stars(paste0("../data_src/oisst/monthly_anom/", i, monj, "_anom.nc")) #Old vers use anom, not sst, so no subtract 30yr mean
      z <- read_stars(paste0("../data_src/oisst/monthly_sst/", i, monj, "_sst.nc"))
      
      names(z)[1] <- "anom" 
      z[[1]][which(landm==1)] <- NA_real_
      ice <- read_stars(paste0("../data_src/oisst/monthly_icemask/", i, monj, "_icemask.nc"))
      z[[1]][which(ice[[1]]==1)] <- NA_real_
      z[[1]] <- z[[1]] - stm[[j]][[1]] #Old-version without this substraction while using _anom directly
      write_stars(z, filet)
    }
  }
}
}

## Just check if plot 2022-04
#if (Recalculate_anom) {
gx2 <- gplotx(z, "anom", minz=-3.5, maxz=3.5, return=TRUE) #This anomaly is substrat long-term mean (30yrs, 1982-2011)
#y <- read_stars("../data_src/oisst/monthly_anom/202204_anom.nc")
#names(y)[1] <- "anom" 
#gy <- gplotx(y, "anom", minz=-3.5, maxz=3.5, return=TRUE)

layt <- rbind(c(1,1,2,2))
grid.arrange(gx2, gy, layout_matrix=layt)

#}

## It's hard/slow to detrend each x,y along almost 30yr (30x12) year trend. Must read 1440x720x30x12 at once...
#removeLargeObj <- TRUE
#if (removeLargeObj) {
#  rm(land_mask, landm, stm)
#}
# 20220727 modified for detrend not distinguish month, all time-span should involved
Reload_allnc <- FALSE
#predyr <- 2021 # previous year that has whole 1-12 month data
#predrng<- seq(1982, predyr) 
yrs <- yrng[length(yrng)] - yrng[1] + 1
#yrs <- predrng[length(predrng)] - predrng[1] + 1

if (Reload_allnc) {
  for (i in yrng) {
    for (j in 1:12) {
      if (i==curryr & j>(currmo-1)) break
      monj <- fifelse(j<10, paste0("0",j), paste0(j))
      jstr <- monstr[j]
      
      z <- read_stars(paste0("../data_src/oisst/monthly_anom_icemask/", i, monj, "_anom.nc"))
      names(z) <- "anom"

      if (i==yrng[1] & j==1) {
        stdrx <- z
        datex<- as.Date(paste0(i, monj, "01"), format="%Y%m%d")
      } else {
        stdrx <- c(stdrx, z)
        datex<- c(datex, as.Date(paste0(i, monj, "01"), format="%Y%m%d"))
      }
    }
  }
  print(paste0("Now in Year-month: ", i," - ", monj, " and have length stdrx: ", length(datex)))
#  save(stdrx, file="../data_src/stats/monthly_anom_1982_2021_nc.RData") #too big, not save
} #else {
#  load("../data_src/stats/monthly_anom_1982_2021_nc.RData")
#}

Recalculate_anom <- FALSE
if (Reload_allnc) {
  trd <- merge(stdrx) %>% st_set_dimensions(3, values = as.POSIXct(datex), names = "time") 
  tx <- st_get_dimension_values(trd, 3)
  
  predtx <- function(x) {
    if (all(is.na(x))) return(list(as.numeric(x)))
    x[!is.na(x)] <- as.numeric(stats::predict(lm(x ~ tx))) 
    return(list(as.numeric(x))) 
  }  
  trdx <- st_apply(adrop(trd), c(1,2), predtx)
  tk <- array(unlist(trdx[[1]]), dim = c(length(datex),1440,720)) 
  save(tk, file="../data_src/stats/predict_anom_1982_202204_nc.RData") #too big, not save
} else {
  if (Recalculate_anom) load("../data_src/stats/predict_anom_1982_202204_nc.RData")
}

if (Recalculate_anom) {
for (k in seq_along(yrng)) {
for (j in 1:12) {
  if (yrng[k]==curryr & j>(currmo-1)) break
  monj <- fifelse(j<10, paste0("0",j), paste0(j))
  jstr <- monstr[j]
  nt <- (k-1)*12+j
#  for (i in yrng) {
#    if (i==curryr & j>(currmo-1)) break
#    
#    z <- read_stars(paste0("../data_src/oisst/monthly_anom_icemask/", i, monj, "_anom.nc"))
#
#    if (i==yrng[1]) {
#      stdrx <- z
#      datex<- as.Date(paste0(i, monj, "01"), format="%Y%m%d")
#    } else {
#      stdrx <- c(stdrx, z)
#      datex<- c(datex, as.Date(paste0(i, monj, "01"), format="%Y%m%d"))
#    }
#  }
#  names(stdrx) <- rep(jstr, length(names(stdrx)))
#  print(paste0("Now in Year-month: ", i," - ", monj, " and have length stdrx: ", length(datex)))
  
#  predy <- function(x) { return (stats::predict(lm(x ~ seq_along(datex)))) }
#  trd <- merge(stdrx) %>% st_set_dimensions(3, values = as.POSIXct(datex), names = "time") %>% 
#    aggregate(by=paste0("months"), FUN = function(y) {
#       if (all(is.na(y))) return(list(as.numeric(y)))
#       y[!is.na(y)] <- as.numeric(predy(y)) #Seems if return value has equal length, the result will be simplified by aggregate if using as.matrix
#       return(list(as.numeric(y))) #Old-version bug: not predy(y), should be y. See previous testing chunk
#    })
#  tk <- array(unlist(trd[[1]]), dim = c(length(datex),1440,720)) #bug fix: NOT each predict have the same length result (if input has NA)
  
#  print(paste0("Predict ok with dim(tk): ", paste0(dim(tk), collapse=",")))
#  for (k in seq_along(yrng)) {
    
    z <- read_stars(paste0("../data_src/oisst/monthly_anom_icemask/", yrng[k], monj, "_anom.nc"))
    #notna1 <- which(!is.na(z[1][[1]]))
    #notna2 <- which(!is.na(tk[k,,]))
    #notna <- intersect(notna1, notna2)
    #z[[1]][notna] <- z[[1]][notna]- tk[1,,][notna]
    z[[1]] <- z[[1]]- tk[nt,,] #Any of z[[1]] or tk[1,,] NA will cause NA here, it makes sense.
    names(z)[1] <- "anom"
    
    filet <- paste0("../data_src/oisst/monthly_anom_detrend/", yrng[k], monj, "_anom.nc")
    write_stars(z, filet)
    print(paste0("Now write Year-month: ", yrng[k], " - ", monj, " and use nth of tk: ", nt, " ok"))
  }
}
}  

yk <- read_stars("../data_src/oisst/monthly_anom/202006_anom.nc")
zk <- read_stars(paste0("../data_src/oisst/monthly_anom_detrend/", 2020, "06", "_anom.nc"))
z <- read_stars(paste0("../data_src/oisst/monthly_anom_icemask/", 2020, "06", "_anom.nc"))

names(yk)[1] <- "anom" 
names(zk)[1] <- "anom"
names(z)[1] <- "anom"
print(range(na.omit(as.data.table(yk)$anom))) #-5.911333 10.598333
print(range(na.omit(as.data.table(zk)$anom))) #-5.39709 16.39870   #old -8.227043 11.727673
print(range(na.omit(as.data.table(z)$anom)))  #-7.724047 10.483911 #old -8.361936  6.181935

##gyk<- gplotx(yk,"anom", minz = -6, maxz = 11.0, returnx=TRUE)
##gz <- gplotx(z, "anom", minz = -6, maxz = 16.5, returnx=TRUE)  
##gzk<- gplotx(zk,"anom", minz = -8, maxz = 11.0, returnx=TRUE)  
#gyk<- gplotx(yk,"anom", minz = -6, maxz = 11.0, returnx=TRUE)
#gzk<- gplotx(zk,"anom", minz = -10, maxz = 7.0, returnx=TRUE)  
#gz <- gplotx(z, "anom", minz = -8, maxz = 11.0, returnx=TRUE)  
gyk <- plot_anom("../data_src/oisst/monthly_anom/202006_anom.nc", returnx=TRUE)
gzk <- plot_anom(paste0("../data_src/oisst/monthly_anom_detrend/", 2020, "06", "_anom.nc"), returnx=TRUE)
gz <- plot_anom(paste0("../data_src/oisst/monthly_anom_icemask/", 2020, "06", "_anom.nc"), returnx=TRUE)

layt <- rbind(c(1,1),c(2,2),c(3,3))
grid.arrange(gyk, gzk, gz, layout_matrix=layt)

#gyk
#gzk
#gz

ga <- plot_anom(paste0("../data_src/oisst/monthly_anom_detrend/", 2022, "04", "_anom.nc"), #zk,
                "anom", #minz = -8, maxz = 11.0, 
                returnx=TRUE) 
ga

if (!require("layer")) {
  if (!require("remotes")) install.packages("remotes"); library(remotes)
  remotes::install_github("marcosci/layer")
  library(layer)
}

Reload_sst_nc <- FALSE
if (Reload_sst_nc) {
  sst1 <- read_stars(paste0("../data_src/oisst/monthly_sst/202002_sst.nc"))
  names(sst1)[1] <- "sst"

  sst202002 <-  as.data.table(sst1) %>% 
    .[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>% 
    #.[,.(longitude, latitude, heatwave)] %>%
    .[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
    .[,`:=`(x=NULL, y=NULL)] %>%
    st_as_sf(coords=c("longitude", "latitude"))

  sst2 <- read_stars(paste0("../data_src/oisst/monthly_sst/201002_sst.nc"))
  names(sst2)[1] <- "sst"

  sst201002 <-  as.data.table(sst2) %>% 
    .[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>% 
    .[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
    .[,`:=`(x=NULL, y=NULL)] %>%
    st_as_sf(coords=c("longitude", "latitude"))

  sst3 <- read_stars(paste0("../data_src/oisst/monthly_sst/200002_sst.nc"))
  names(sst3)[1] <- "sst"

  sst200002 <-  as.data.table(sst3) %>% 
    .[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>% 
    .[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
    .[,`:=`(x=NULL, y=NULL)] %>%
    st_as_sf(coords=c("longitude", "latitude"))

  sst4 <- read_stars(paste0("../data_src/oisst/monthly_sst/199002_sst.nc"))
  names(sst4)[1] <- "sst"

  sst199002 <-  as.data.table(sst4) %>% 
    .[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>% 
    .[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
    .[,`:=`(x=NULL, y=NULL)] %>%
    st_as_sf(coords=c("longitude", "latitude"))
  
  save(sst199002, sst200002, sst201002, sst202002, file="../data_src/stats/sst_1990_2020_stack.RData")
} else {
  load("../data_src/stats/sst_1990_2020_stack.RData")
}

tl0 <- layer::tilt_map(sst199002)
tl1 <- layer::tilt_map(sst200002, y_shift=75)
tl2 <- layer::tilt_map(sst201002, y_shift=150)
tl3 <- layer::tilt_map(sst202002, y_shift=225)

map_list <- list(tl0, tl1, tl2, tl3)

#palettem: https://cran.r-project.org/web/packages/viridis/vignettes/intro-to-viridis.html
layer::plot_tiltedmaps(map_list, palette = "turbo")

zt1 <- fread("../data_src/verify/198201_detrend_anom_yeh.csv") 
zt <- read_stars(paste0("../data_src/oisst/monthly_anom_detrend/", 1982, "01", "_anom.nc"))
tt1 <- as.matrix(zt1)
tt <- zt[[1]]
quantile(na.omit(tt1[25,]))

##          0%         25%         50%         75%        100% 
## -1.57301116 -0.26347324  0.04314169  0.33155985  1.47273052

quantile(na.omit(tt[25,]))

##           0%          25%          50%          75%         100% 
## -1.866937160 -0.298840664  0.005112544  0.310795255  1.452879071

length(which(is.na(tt1)))

## [1] 548586

#[1] 548586
length(which(is.na(tt))) #[1] 474337 !! large diff, should check...

## [1] 474337

range(na.omit(as.numeric(tt1)))

## [1] -7.042250  3.954143

range(na.omit(as.numeric(tt)))

## [1] -7.016757  3.966041

zt <- as.data.table(zt) %>% setnames(3, "anom")
xval <- sort(unique(zt$x))
yval <- rev(sort(unique(zt$y)))
zt1 <- cbind(c(1:1440), zt1)
setnames(zt1, 1, "xid")
zt1a <- melt(zt1, id.vars = "xid", measure.vars = patterns("^V")) %>%
  .[,vart:=gsub("V","",variable) %>% as.integer()] %>%
  .[,`:=`(x=xval[xid], y=rev(yval)[vart])] %>%
  .[,.(x, y, value)] %>% setnames(3, "anom")

setorder(zt, -y, x) ## the same order as original zt data.frame
setorder(zt1a, -y, x)

gz1a <- gplotx(zt1a, "anom", minz = -7.5, maxz = 4, returnx=TRUE)
gzt <- gplotx(zt, "anom", minz = -7.5, maxz = 4, returnx=TRUE)
layt <- rbind(c(1,1),c(2,2))
grid.arrange(gzt, gz1a, layout_matrix=layt)

plot(1:720, zt1a[x==121.125,]$anom, type="l", col="blue")  #almost the same
points(1:720, zt[x==121.125,]$anom, type="l", col="red", lty=2, add=T)

plot(1:720, zt1a[x==10.125,]$anom, type="l", col="green")
points(1:720, zt[x==10.125,]$anom, type="l", col="brown", lty=2, add=T)

chk_x <- 360-145+0.125 ## Now try to plot paper 145W, 55N detrend vs anom not detrended
chk_y <- 50.125
plot(1:720, zt1a[x==chk_x,]$anom, type="l", col="black")
points(1:720, zt[x==chk_x,]$anom, type="l", col="magenta", lty=2, add=T)

plot(1:1440, zt1a[y==chk_y,]$anom, type="l", col="black")
points(1:1440, zt[y==chk_y,]$anom, type="l", col="magenta", lty=2, add=T)

plot(1:1440, zt1a[y==0.125,]$anom, type="l", col="black")
points(1:1440, zt[y==0.125,]$anom, type="l", col="magenta", lty=2, add=T)

Reverify_paper <- FALSE
idxt <- which(zt$x==chk_x & zt$y==chk_y) #229821
zt[idxt,] #-0.3095802

##          x      y       anom
## 1: 215.125 50.125 -0.3095802

idxm <- idxt %% 1440
idym <- floor(idxt/1440) +1
tt[idxm,idym] == zt[idxt,]$anom 

## [1] TRUE

if (Reverify_paper) {
  anomx <- c()
  anomd <- c()
  sstx <- c()
  for (i in yrng) {
    for (j in 1:12) {
      if (i==curryr & j>(currmo-1)) break
      monj <- fifelse(j<10, paste0("0",j), paste0(j))

      zd <- read_stars(paste0("../data_src/oisst/monthly_anom_detrend/", i, monj, "_anom.nc"))
      anomd <- c(anomd, zd[[1]][idxm, idym])
      
      zk <- read_stars(paste0("../data_src/oisst/monthly_anom_icemask/", i, monj, "_anom.nc"))
      anomx <- c(anomx, zk[[1]][idxm, idym])

      zs <- read_stars(paste0("../data_src/oisst/monthly_sst/", i, monj, "_sst.nc"))
      sstx <- c(sstx, zs[[1]][idxm, idym])
      if (i==yrng[1] & j==1) {
        testdatex <- as.Date(paste0(i, monj, "01"), format="%Y%m%d")
      } else {
        testdatex<- c(testdatex, as.Date(paste0(i, monj, "01"), format="%Y%m%d"))
      }
    }
  }  
  save(anomx, anomd, sstx, testdatex, file="../data_src/stats/verify_paper_loc145W50N_sst_anom_1982_202204.RData")
} else {
  load("../data_src/stats/verify_paper_loc145W50N_sst_anom_1982_202204.RData")
}

plot(testdatex[1:468], sstx[1:468], type="l", col="magenta")

plot(testdatex[1:468], anomx[1:468], type="l", col="magenta", ylim=c(-2,3))
points(testdatex[1:468], anomd[1:468], type="l", col="black", add=T)